function [J, dxi] = rcnl_jacobian(ddelta, x, z, W, fe, z_alt)

    nfe    = size(fe, 2);
    xx     = fe' * fe;
    x1     = x(:, 1:end-nfe);
    xtilda = x1 - fe * (xx \ (fe' * x1));

    XZ       = xtilda' * z_alt;
    A        = inv(XZ * W * XZ') * XZ * W * z_alt';
    J        = ddelta' * (z' - (z' * xtilda) * A)';
    dxi      = (ddelta' - (ddelta' * A') * xtilda')';
    dxi      = dxi - fe * (xx \ (fe' * dxi));
end
